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We present developing of method of the numerical analysis of polarization in the Gauss- 
Legendre Sky Pixelization (GLESP) scheme for the CMB maps. This incorporation of 
the polarization transforms in the pixelization scheme GLESP completes the creation of 
our new method for the numerical analysis of CMB maps. The comparison of GLESP 
and HEALPix calculations is done. 

1. Introduction 

The analysis of the anisotropy of CMB temperature and polarization is one of the 
most effective methods for the extraction of the cosmological information and for 
tests of cosmology and fundamental physics. 

The measurements of the anisotropy of CMB temperature performed by the 
WMAP mission allow us to establish parameters of the cosmological model of the 
Universe with unprecedented precision Refs. 1, 2, 3, 4. 

Increasing sensitivity and angular resolution of the CMB data, including recently 
available WMAP 5 year data Ref. 2 , ACBAR Ref. 5 QUaD Refs. 6,7,stimulates 
significant development and increasing predictability of the corresponding software 
(see, for instance, CAMB Ref. 8, COSMOMC Ref. 9, RICO Ref. 10, etc.) 

During the next decade, after the PLANCK experiment, an investigation of the 
CMB polarization (including the B-mode) will be at the focus of the CMB science. 
Planning the CMBpol mission Ref. 11, B-pol mission Ref. 12, etc. requires significant 
improvement in estimation of the errors of the signal. Partially an uncertainties due 
to pixelization of the CMB sky and their propogation to the CMB power and the 
map can be potential sources of the error. For planning high resolution the CMB 
experiments the most frequently used the HEALPix package Ref. 13 as well as 
GLESP Ref. 14, ECP Ref. 15, and some others, needs to be tested in order to 
provide an exact information about the error bars of the convolution of the T^E^B- 
maps to the corresponding coefficients of the spherical harmonics decomposition not 
only for the power spectrum, but for the real and imaginary part of the coefficients 
as well. 

The information about the multipole structure of the CMB signal is vital for the 
low multipoles (^ = 2, 3..), since a lot of theoretical predictions about the properties 
of the cosmological model are related directly to the global morphology of the 
signal. As an example, we would like to mention widely discussed the Bianchi Vllh 
anisotropic cosmological model Refs. 16, 17, 18 which can mimic the anisotropy of 
the CMB power at the range of multipoles ^ < 20, and the Cold Spot Ref. 19 as 
well. The whole sky decomposition is vital for testing the alignment and planarity 
of the CMB anisotropy multipoles at 2 < ^ < 5, discussed in Ref. 20 The exact 
information about the phases of the CMB signal is very useful for investigation of 
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statistical anisotropy and non-Gaussianity of the CMB Refs. 21, 22, 23, 24. 

This Paper is devoted to presentation of a new package GLESP-pol (see Ap- 
pendix), and investigation in details the errors of the standard transition "map to 
a^^^" and vice verse for the most frequently used HEALPix2.11 and newly released 
the GLESP-pol package. We would like to point out that both these packages re- 
veal some peculiarities of the reconstruction of the coefficients of decomposition, 
especially for polarization. The major part of the error belongs to the ^, m = 0, 
and £,m = 2 modes, when the simplest variants of decomposition were used. For 
the HEALPix2.11, it is "zero iteration" key which blocks the correction of the ai^m 
taken from the map. For the GLESP-pol package, the maxima of the ai^rn errors 
correspond to the GLESP 1.0 pixelization Refs. 14, 25. However, all these problems 
can be successfully resolved by implementation of iterations for the HEALPix2.11 
(the key "iterative analysis " 3 or 4 iterations) and the GLESP-pol pixelization for 
the polarization. 

The outline of the paper is the follows. In Section 2, we discuss the difference 
between the HEALPix and the GLESP-pol scheme of the CMB sky pixelization, 
focusing on the polarization of the CMB. Section 3 is devoted to investigation of the 
errors of the "map^ ai^m map" transition for the CMB temperature anisotropy 
for the HEALPix and the GLESP-pol. In Section 4, we discuss the same issue for 
the Q,U Stokes parameters and E and B-modes of polarization. In Appendix the 
basic relations and description of the GLESP-pol package are presented. 

2. Basic definitions 

The temperature and polarization CMB anisotropy can be described in terms of the 
Stokes parameters T^Q^U through spherical harmonics decomposition Yi^^iO^cj))^ 
and spin ±2 spherical harmonics±2^,m(^5 0) 

^ m 

Q{0, cj)) ± iU{e, <t>) = Y.Y. ±2a^,m±2>£,m(^, <l>) (1) 

I m 

Here ^, (j) are the polar and azimuthal angles of the polar system of coordinates, 
ai^rn stands for the temperature anisotropy and the spin coefficients ±2<^^,m can be 
decomposed into E and B modes of polarization (see Appendix for details): 

±2a,,m = -{aim ± i4m), <f = {-^T {^^ )l-m (2) 

The conversion of the T^Q^U signals to corresponding ai^rn and ±2Cii,m coeffi- 
cients is given by following integrals: 

j dx j dct> (r(x, 0), Q(x, 0), U{x, 0)) o,±2F,:^(x, 0), 

(3) 
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where index marks the temperature anisotropy, and x = cos^. As it seen from 
Eq(3), the mathematical basis of any schemes of the pixelization of the CMB sky 
is very simple. We need to estimate the integrals in Eq(3) with very high accu- 
racy taking into account the properties of discrete representation of the signal on 
the sphere. However, the modern CMB experiments normally deal with the incom- 
plete sky due to peculiarities of design, scan strategy or implementation of different 
sort of mask. In this case the scheme of the pixelization of the sky become even 
more important, since we need to use the pixel domain for estimation of the power 
spectrum and investigation of the statistical properties of the CMB signal without 
implementation of the o,±2^i^,m-coefficients. For these purpose the basic idea of the 
HEALPix package (equal area isolatitude pixelization) is very useful and more ad- 
vanced in comparison to other pixelization. However, for the whole sky analysis of 
the T, Q, U signals o,±2^i^,m-domain seems to be more optimal, in terms of the CPU 
timing, as from the scientific point of view. This is why in this paper we propose 
the GLESP-pol pixelization as complementary approach to the HEALPix package. 

For numerical evaluation of the integral Eq(3), we use the Gaussian quadratures 
Ref. 14. This approach was proposed by Gauss in 1814, and developed later by 
Christoffel in 1877. As the integral over x in Eq. (3) is an integral over a polynomial 
of X we can use the following equality Ref. 26: 

/ dx^T{x, mm{^. 0) = 5^^iAr(x,-, 0) . (4) 

where both AT{xj ^ (j))Y/^{xj , (/)) and the proper Gaussian quadrature weighting 
functions, Wj = w{xj), are taken at points xj which are the net of roots of the 
Legendre polynomial 

Pn{xj) = 0. (5) 

Here N is the maximal rank of the polynomial under consideration. 

It is well known that the equation P/v {xj ) = has number of zeros in interval 
— l<x<l. For the Gaussian-Legendre method Eq(4), the weighting coefficients 
are 

2 / 

where ' denotes a derivative. They can be calculated together with the set of xj 
with the 'gauleg' code Ref. 26. 

The old GLESP (version 1.0) pixelization scheme Ref. 25 was defined as follows: 

• In the polar direction x = cos 0, we define Xj, j = 1, 2, . . . , A/", as the net of 
roots of Eq. (5). 

• Each root Xj determines the position of a ring with A^^ pixel centers with 
^-coordinates 0^. 

• All the pixels have nearly equal area. 

• Each pixel has weight wj (see Eq (6)). 
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This scheme for the teraperature anisotropy was reahzed in the pubhcly available 
code presented in www . glesp . nbi . dk and below we call it "GLESP 1.0" denoting as 
'grA' (grid of equal areas). The new code GLESP-pol is based on the same definition 
of the roots of the Legendre polynomials, but different definition of the pixel area 
(see the item 3 from the top): 

• All the rings have the same number of pixels (case 'grN' — grid of equal 
number of pixels in ring) by default. That is a cylindric projection of the 
sphere. We have also checked a special case when the number pixels of ring 
has an increment 4 starting from 10 pixels near the poles (case 'grS' — 
special grid). The old scheme (GLESP 1.0) when all the pixels have nearly 
equal area (grA) is also accessible. 

In fig.l, we show the diff'erences in pixelization of the GLESP 1.0, the GLESP- 
pol (grN), the GLESP-pol (grS) and the HEALPix. One can see, in the central 
part of the map the properties of the GLESP-pol (grN,grS) are the same as for the 
GLESP 1.0, but in the vicinity of the polar cups they are significantly different due to 
over-pixelization. This modification is vital for estimation of the o ±2^^^ m~ 

coefficients 

due to specific behaviour of ±2^^,771(^5 0) -spherical harmonics (see Appendix for 
details). 

3. Errors "ai^rn — ^ map — > CLi^rn'' transition 

This section is devoted to estimation of the error of the ai^rn reconstruction by 
implementation different regimes of the HEALPix 2. 11 and the GLESP-pol (grN 
and grS) packages. The error bars were defined in the following way. Let us take the 
coefficients aj'^'^ for temperature anisotropy T and E and B modes of polarization 
for ACDM concordance model by implementation of the Monte Carlo simulation of 
the random Gaussian signal^. For the C^s-spectrum, we have taken a value equal 
to 10~^^ just to escape zero in the test calculations. 

Then, by implementation of the HEALPix 2. 11 and the GLESP-pol packages 
we used these ^ ' coefficients to create the map of the signal, keeping for both 
packages the same number of pixels. Thus, the HEALPix map M|^, and the GLESP- 
pol map Mq are defined as 

Mi,=Hal'^'^, M^ = Gal'^'^ (7) 

where H and G are the HEALPix and the GLESP-pol operators for "a^,r7i to map" 
transition, the index i marks the number of iteration (0 as default, or 1-4), the index 
p = 1 marks the GLESP-pol grN pixelization, and p = 2 corresponds to the grS 
pixelization. Let us define the corresponding transition "map to a^,r7i" as for 



^We have used the Gaussian signal for simpUcity. However, only for the Gaussian signal the power 
spectrum C{£) = {2£ + Em=-£ I^Jm^l^ unique characteristic, which determines all 

the statistical properties. 
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Fig. 1. Left column is the Molldweide projection of pixelization grids (from top to bottom): (a) 
the standard pixeUzation grid of the GLESP 1.0 (b) the GLESP-pol rectangular grid with the same 
number of pixels per each ring (the so called case 'grN'), (c) the GLESP-pol grid with an pixel 
number increment 4 (case 'grS') starting from 10 pixels near poles but not greater than a given 
resolution in the equator ring, (d) the HEALPix grid. The right column shows the corresponding 
pixelization in the vicinity of the polar cups. 

the HEALPix and G"^ for the GLESP-pol: 

<f = G-Mg = G-Ga,^f (8) 

where bj'^'^ and cj^'^ are now the reconstructed coefficients for the HEALPix 
and GLESP-pol correspondingly. For idealistic case, when G£ ^ — > map — > ag jj^'' 

T E B T E B 

transition has no error bars, the reconstructed 6^ ^ ' and ^ ' have to be identical 
to the input coefficients o^f and R-^H = HH"^ = I, G'^G = GG"^ = I, 
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where I is just a unit matrix. In reality, neither the HEALPix, nor the GLESP- 
pol packages have non-zero error of the reconstruction due to window functions 
of the pixels and computational errors for the spherical harmonics. That means 
that corresponding absolute errors for the " ai^rn map^ tt^,m" transition can be 
defined as follows: 

(9) 

where 5Re and stand for the real and imaginary parts of the coefficients. Thus, 
the relative error is given by 

pi ri 
^^Km ) ^^(^Am ) 

DP JP 

p G£,m p Gi,m 

^^y^l^m ) ^^K,m ) 



(10) 



Note that defined in Eq(lO) relative errors are related to the error of the power 
spectrum 



2 



where gi^rn and a,^^ denote the reconstructed and input multipole coefficients, 
and g£^rn = ci£,m{^ + ^£,m)- Taking into account that 5£^rn ~^ ^1 m — 25Re(^,,^) and 
^e{Si^rn) = [{^eai^m)'^ri^rn + {^ma£^rnYy£ 

,m]/|^^,m| 5 where 777, and yi rn denote 
relative errors for real and imaginary parts from Eq(lO), we gets: 

C£ ~ EmM' ^ ' 

From Eq(12) clearly seen that the error of reconstruction of real and imaginary 
part of each ^, m coefficient propagates to the error of the power spectrum through 
weighting coefficients 

(Jiea,,^)^ ^ y {^ma£ 



<m = ^ r- . (13) 



and formally depends on the power spectrum and the morphology of the input 
signal. This is why in addition to the input random Gaussian CMB signal we will 
discuss the errors of reconstruction of the multipoles and the power spectrum for 
very asymmetric maps, like WMAP 5 Q and V bands. 
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3.1. Errors for HEALPix and GLESP-pol temperature anisotropy 

We will start our analysis from estimation of the dynamical range of variations of 
the temperature anisotropy for Nside=1024. For that we will use the input signal, 
which corresponds to random Gaussian CMB, and Internal Linear Combination 
Map (WILC5) from the LAMBDA archive (http://lambda.gsfc.nasa.gov/) and 
corresponding K, Ka, Q, V and W total channels maps. In Fig2 we plot the diagram 
A{i) versus A{th), whereA(z) = string {\a£^m)\^ where the operator "string" trans- 
forms the \a£^m\ coefficients to one dimensional string |ai,oU<^i,i|--- for K-W bands, 
and A(th) is the string for the random Gaussian CMB signal. From this diagram 




10-' 10"'' IQ-" 10-' 10-* 10-' 10-^ 10"' 10-8 ^0-5 

A(th) A(th) 1 



Fig. 2. Left. The A{i) versus A{th) diagram for the WILC5 (the black dots), K-band (the blue 
dots), KA-band (the green dots), V-band (the red dots). Both A(i) and A(th) are in mK. Middle. 
The same as the left plot, but for replacement of the WILC5 by A(th). Right. The power spectrum 
C{i)mK'^ for the WILC5 (the red bottom line), for random realization (the black line), the K band 
( the solid blue line), the KA-band (the green solid line), and for V band (the red solid line). 

one can see that, for example, to estimate the a£^rn coefficients by using ILC method 
we should have, at least, the relative error better than 10~^ — 10~^ if the low fre- 
quency K-band is included to the analysis. T he accuracy of the a^^rn coefficients 
should be about 3-4 orders of magnitude better if we are interested in different sort 
of coupling between different multipoles (so called the non-Gaussianity tool). 

To estimate the errors of the ai^m map^ <^£,m" transition, let us start from 
the analysis of the temperature anisotropy for the HEALPix 2. 11 and the GLESP- 
pol packages. In Fig. 3, we show the maps for ^ — a^^^-signal, where ^ obtained 
by implementation of the HEALPix 2. 11 "zero iteration" key (top left map) and 4 
iterations (top right map). 

One can see, that the top left map reveals all the peculiarities of the pixelization, 
localized at the vicinity of the North and the South poles. The top left plot show the 
map of differences, but after 4 iterations. No visible large scale defects can be found. 
The bottom left and right maps represent the GLESP-pol pixehzation without any 
iterations for the same number of pixels, as the top one. In Fig. 4, we plot the real 
and imaginary parts of the absolute errors from Eq(9) for the maps, shown in fig. 3). 

The definition of function A{in) and r in Fig. 4 is the following: A{in) = 
string {a£^rri)^ where the operator "string" transforms the a£^m coefficients to one 
dimensional string ai,o, and r = k?^^,^, ^If^,^, (see Eq.(lO). 

Fig. 4 clearly show that the HEALPix 2. 11 (no iteration) reconstruction is differ- 



April 16, 2009 19:19 WSPC/INSTRUCTION FILE glesp-ijmpd 




Fig. 3. The differences of reconstructed and input maps for the HEALPix2.11 (the top pair) 
and the GLESP-pol (the bottom pair). Top left plot corresponds to "zero iteration" key (the 
color scale is — 10~^, 10~^mK), top right is for 4 iterations (the color scale is — 10~^, 10~^mK). 
Bottom left-the GLESP-pol reconstruction for the grN mode (the color scale is — 10~^, 10~^mK). 
Top right corresponds to the grS mode (the color scale is — 10~^, 10~^mK). The number of pixels 
for the HEALPix2.11 and the GLESP-pol are practically the same. No correction by the window 
function of the pixels. 





Fig. 4. The errors of reconstructions for the HEALPix2.11 (the top pair) and the GLESP-pol 
(the bottom pair). Top left plot corresponds to "zero iteration" key, top right is for 4 iterations. 
Bottom left — the GLESP-pol reconstruction for the grN mode. Top right corresponds to the grS 
mode. Black dots corresponds to the real part of ^, m modes, the red dots show for the imaginary 
part. 



ent for the real and imaginary parts of the signal. For real part (see Fig. 4, top left), 
the relative error has the secondary zone, localized at r > 10~^, while for imagi- 
nary part the major part of the points belongs to r < 10~^. Note that even after 
implementation of 4 iterations for the HEALPix2.11 and grN, grS modes of the 
GLESP-pol, there exists a very small number (~ 10) of modes with error 1%-10%. 
Thus, the main conclusion is that for Nside=1024 the HEALPix iterations signifi- 
cantly improve the global morphology of the map and the error of reconstruction. 
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which is practicahy the same the HEALPix and the GLESP-pol. 

The next question, which we would like to discuss, is how the error of recon- 
struction depends on the number of pixels and their size. Both these parameters are 
determined by the choice of Nside for the HEALPix, or by the maximal resolution 
of the map imax for the GLESP-pol. Answering this question in Fig. 5 and Fig. 6, 
we plot the corresponding maps of errors and A{in) versus r diagrams, similar to 
Fig. 3 and Fig. 4, but for less number of pixels. 




Fig. 5. The differences of reconstructed and input maps for the HEALPix 2.11 for Nside=32 (the 
top pair) and the GLESP-pol (the bottom pair). Top left plot corresponds to "zero iteration" key 
(the color scale is — 10~^, 10~^mK), top right is for 4 iterations (the color scale is —8 • 10~^,8 • 
10~^mK). Bottom left — the GLESP-pol reconstruction for the grN mode (the color scale is 
— 10~^, 10~^mK). Bottom right corresponds to the grS mode (the color scale is —2.5 • 10~^, 2.5 • 
10~^ mK). The number of pixels for the HEALPix 2.11 and the GLESP-pol are practically the 
same. No correction by the window function of the pixels. 



It would be important to note that as for the high resolution map of difference, 
shown in Fig. 3, as for the low resolution map (see Fig. 5) for the HEALPix 2.11 with 
zero iteration the major component of the error is related to 6^,rn=o mode. These 
harmonics manifest themselves as horizontal lines parallel to the Galactic plane. 
Two horizontal lines with high amplitude signal along mark the HEALPix zones, 
where the number of pixels for each equal latitude ring start to decrease, when 6^0 
(the North pole), or ^ ^ tt (the South pole). However, after 4 iterations all these 
peculiarities of the map of errors were significantly suppressed, except two zones 
around the North and the South cups. For the GLESP-pol package with imax = 32 
the minimal level of errors for the reconstructed map is given by the grN pixelization, 
when r <C 10~^ for major part of the pixels, and r <C 10~^ for the grS pixelization. It 
would be important to note that the accuracy of reconstruction of real and imaginary 
parts of the coefficients of expansion are different for the HEALPix and the GLESP- 
pol. As it is follows from Fig4 and Fig6, the HEALPix recovers the imaginary part 
significantly better than a real one even for zero iteration mode. The GLESP-pol 
reconstructs real and imaginary parts with nearly equal errors. 
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0.004 - - 0.004 - 

0.002 - . - 0.002 - 

0.000 - r:'" '• - ^ 0.000 - 

-0.002 - ' - -0.002 - 



Fig. 6. The errors of reconstructions for the HEALPix2.11 (Nside=32) (the top pair) and the 
GLESP-pol (the bottom pair). Top left plot corresponds to "zero iteration" key, top right is for 4 
iterations. Bottom left — the GLESP-pol reconstruction for the grN mode. Top right corresponds 
to the grS mode. Black dots corresponds to the real part of i,m modes, the red dots are for 
imaginary part. 



4. The GLESP-pol and the HEALPix polarization. 

As we have mentioned already in Introduction, the GLESP-pol pixelization was 
designed to assess the problem of accurate reconstruction of the coefficients of spin 
±2 spherical harmonics decomposition (see Appendix for details). Since the spin 
±2 spherical harmonics have a peculiarity, from the computational point of view, 
behaviour in the vicinity of the polar cups, the differences in the GLESP-pol po- 
larization and the HEALPix produce different error of reconstruction and become 
more visible especially for polarization. In Fig. 7, we plot the map of differences 
between input and output signal, reconstructed by the GLESP-pol (grN and grS 
pixelization). 

As one can see from this figure, the reconstruction of Q and U components by 
the grN and grS pixelization is characterized by very high accuracy, but the grN 
pixelization looks slightly better. For the HEALPix 2.11 the corresponding maps for 
differences are shown in Fig. 8. 

For low resolution pixelization with Nside=32 and imax = 32 the corresponding 
maps are shown in Fig. 9 

Thus, one can see that implementation of the 4 iterations for the HEALPix 
package gives us practically the same result, as from GLESP-pol even for the low 
resolution maps (see Fig. 10). In Fig. 11 we show the diagram, similar to the Fig. 4, 
but for E and B components of polarization. 

One can see that the GLESP-pol mode grN gives us practically the same errors, 
as the HEALPix 2.11 4 iteration mode. However, the CPU time for the GLESP-pol 
is about 4 time smaller then for the HEALPix 2.11 due to absence of iterations. 
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Fig. 7. The differences of reconstructed and input maps for the GLESP-pol for £max = 1500 (the 
top pair is for Q the Stokes parameter) and the bottom pair is for U). Top left plot corresponds to 
the grN pixelization (the color scale is —3 - 10"-'^°, 3- 10~^^mK), top right is for the grS pixelization 
( the color scale is — 10~^, 10~^mK). Bottom left — the GLESP-pol reconstruction of Q for grN 
mode (the color scale is —2.5 • 10"^*-*, 2.5 • 10~^^mK). Bottom right corresponds to the grS mode 
(the color scale is — 10~^, 10~^mK). 




Fig. 8. The differences of reconstructed and input maps for the HEALPix2.11 for Nside=1024 
(the top pair corresponds to Q components and the bottom pair is for U. Top left plot corre- 
sponds to the "zero iteration" key (the color scale is — 10~^, 10~^mK), top right is for 4 iterations 
(the color scale is 10"^*-", 10~^*-*mX). Bottom left — the GLESP-pol reconstruction for the grN 
mode (the color scale is 10~^, 10~^mK). Bottom right corresponds to the grS mode (the color 
scale is -10-^°, IQ-^^mK). The number of pixels for the HEALPix2.11 and the GLESP-pol are 
practically the same. No correction by the window function of the pixels. 

Finally, we demonstrate absolute relative accuracy for power spectra 
{AC{£)/C{£)) for the GLESP2.0 (grS and grN) and for the HEALPix2.11 (0 and 
4 iteration modes) shown in Fig. (12). As one could see the realtive accuracy for 
restoration of the power spectrum of temperature anisotropy in the GLESP is ap- 
proximately the same for both types grS and grN, and the HEALPix 4 iteration 
mode give better accuracy at low multipoles and reaches the GLESP one at higher. 
For polarization, we have approximately the same accuracy for the GLESP grN and 
HEALPix 4 iteration modes. 
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Fig. 9. The differences of reconstructed and input maps for the HEALPix2.11 for Nside=32 (the 
top pair corresponds to Q components and the bottom pair is for U. Top left plot corresponds 
to the "zero iteration" key (the color scale is — 10~^, 10~^mK), top right is for 4 iterations (the 
color scale is —3 • 10"^*-', 3 • 10~^^mK). Bottom left — the U component for the zero iteration 
(the color scale is — 10~^, 10~^mK). Bottom right corresponds to 4 iterations (the color scale is 
-3 • 10-10,3 • 10-iOmK). 



Fig. 10. The errors of reconstructions for the GLESP-pol (Q for the top pair and U for the 
bottom pair). Top left plot corresponds to the grN pixelization (the color scale is —2 • 10-^^,2 • 
IQ-^^mK), top right is for grN (the color scale is -10"^, lO-^mK). Bottom left — the GLESP- 
pol reconstruction for the grS mode (the color scale is —2 • 10--*^^, 2 • lO^-'^^mK. Bottom right 
corresponds to the grS mode (the color scale is —10"^, 10~^mK). 



5. Conclusions 

Here we present the GLESP-pol package which incorporates calculations of polar- 
ization on the sphere into the CMB analysis package based on the Gauss-Legendre 
Sky Pixelization. We developed corresponding software for data processing. 

According to our numerical calculations, the described scheme for polarization 
preserves the same level of precision as GLESP 1.0 can provide for the temperature 
anisotropy. 

The code continues to be developed now and different definitions of polarization 
modes Refs. 27, 13 are planned to be included. GLESP-pol is open for new ap- 
proaches and could be implemented for extension of the present code, e.g. such as 
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Fig. 11. The errors of reconstructions of ^eaf^ and ^maf^. Top row corresponds to the 
HEALPix2.11 (Nside=1024). From the left to the right are: E mode with zero iteration, E mode 
with 4 iteration, B mode with zero iteration, B mode with 4 iteration. Bottom row is the same, 
but for the GLESP-pol grN and grS modes and £max = 1500. Black dots corresponds to the real 
part of i, m modes, the red dots are for imaginary part. 




500 1000 1500 500 1000 1500 500 1000 1500 



Fig. 12. The relative accuracy of C{£) specrum restoration for temperature anisotropy (left), 
E-polarization (middle) and B-polarization(right). The grS GLESP type is shown with the black 
line, the grN type is shown with red, the HEALpix calculations of C(i) of iteration mode are 
plotted with green color, and of 4 iteration with blue. 



the fast spin-weighted harmonics calculation Ref. 28. A completely new algorithm 
for fast Spherical Harmonics Transform has been proposed recently by Mark Tygert 
Ref. 29. It is worth mentioning that this algorithm can be applied only for GLESP 
pixelization of the sky and effective as 0(A x £'^Log2i) operations instead of usual 
0{i^). Unfortunately a huge prefactor A makes this algorithm only about 3 times 
faster than existing HEALPix and GLESP for £ = 2048. 

The most important part of our investigation is that the HEALPix 2. 11 package 
provides the same accuracy as the GLESP-pol only by implementation of 4 itera- 
tions. The zero iteration mode of the HEALPix 2.11 can provide significant error 
for the coefficients of expansion and be used with caution. 
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6. Appendix 

6.1. Description of polarization 

6.1.1. Basic relations for the flat space 

The polarization is characterized by the symmetric traceless matrix T composed by 
two Stokes parameters, Q(x) &t/(x): 

U 

The functions Q and U depend upon the coordinate frame and components of the 
tensor T- obey the corresponding tensor transformation law: 

f^=oioin (15) 

where the coordinate transformation is given by = 0\x^ . In particular case under 
rotation of the coordinate system with the rotation angle (j) 

cos,/, sin ^\ 
\ — sm cos J 

the functions Q and U are transformed as: 

Q' = Qcos20 + ?7sin20 
U' = -Q sin 20 + t/ cos 20 (17) 

The intensity P = Q'^ -\-U^ = \T^Tl does not depend upon the coordinate system. 

6.1.2. Basic relations on the sphere 

An arbitrary traceless symmetric tensor can be presented in terms of scalar "elec- 
trical", and pseudo-scalar "magnetic", 5, potentials as: 



(18) 



where eij is the completely antisymmetric tensor (en = €22 = 0, ei2 = —621 = 
—sinO) and ';' denotes covariant differentiation on the two dimensional sphere with 
the metric 

dl'^ = dO^ + sin^ Odif^ 
Here 0, (p are the polar and azimuthal angles. In the case we get instead of (14): 
= -T| = Q T^ = sin(9t/, Tl = U/ sini9 , 

and again we have 

/2 = O.bT^^Tl = Q^^U^. 
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Stokes parameters Q and U are linked with the potentials E and B as follows: 

Q = DiE - D2B, U = D2E + DiB , (22) 
and operators Di and D2 are: 

^ . 1 cote d\ ( X d 

Do = 



sin 6 dOdif sin dcp J \ dxdcp 1 — x'^ dip J 

where x = cosO. These operators are identical with expressions (2.22) and (2.23) in 
Ref. 27. 



6.1.3. Spherical harmonics 

Following to Refs. 14, reficiteglespa we use the representation of scalar ("electrical"), 
E^ and pseudoscalar ("magnetic"), 5, polarization potentials identical to (1): 



E = 



B = 



^ E + 2 E fn^Mm cos(m^) - sin(m^)] , (25) 

£=2 \ m=l / 

-I ^max ( ^ \ 

^ E + 2 E /r(a^)[&?„, cos{m(^) - sin(m(^)] , (26) 



/r = y(^+o.5)|^-^Br, o<m<^ (27) 

Here BJ^{x) and fp{x) are the associated Legendre functions (ordinary and normal- 
ized) and Cirn and are the coefficients of decomposition characterizing properties 
of polarization. 

4m = J ^'^^ Jo '^)frix) cosnvfi, (28) 

^= / (ix / d(pE{x,(p)fY^{x) sin m(p . 
2ti J -I Jo 



Similar expressions for the "magnetic" mode can be obtained by replacing of 
^im & ^im with bj^ & bf^ and E with 5. 

As in the previous package the functions fp{x) are found recursively: 



/r(-) = ^i^fr-, - Vg^ '%.'!;r' /r-, , (30) 



or 

rm — 1 



i'm V / J t / ' ' £m—z /Q-j \ 



2(771-1) xf'l'^' /£ + 2-m £ + m-l „_2 
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2(2m)!! 



m+l 



(32) 



Relation (30) starts with /™ and (32) and generates functions /™ for all 

£ > m. Relation (31) starts with /™ and and generates functions /™ for all 
1 < m < Z. In the case function /° must be found with relation (30). 
It is easy to see that 



DiE = F^lel^ cos{mip) - e|„ sm{mip)] 
D2E = ^T[4m sin(m<p) + e|„ cos(m(^)] , 



(33) 
(34) 



where functions F'p and are expressed through the normahzed Legendre func- 
tions (30,31): 



2£+l 
21-1 



1 — X' 



1-^2 



J£ 5 



2P 1 

^(^2_„2)ym^_(^_l)^J« 



or by other way 



m — 1 
l-a;2 



£2 + £ - 2m 



p, _ 2 _|_ £ _ ^ 



l-a;2 



-m 



v/(£ + m+l)(£-m) I IN a; 



Here 



= 0.25{e + 2){£+ !)£{£ - 1), 



(35) 
(36) 

(37) 
(38) 

(39) 



and functions Fp and are normahzed by condition 



In particular. 



m — 1 1 -\- X 
-m 



2Mm "^l-x^*^^' 
m(m - 1) 



rnxfl^ 



, V2m + 3 



m+l 



m — 1 m + l 



l-x2 



Combining (22), (23), (25) and (26), we get for functions Q, & t/: 



/7 
n 



. V2m + 3 1 - mx^ 



(42) 



-= / d^Q cos = ^ (Fr(^)e?™ - *r(^)&lm) , (43) 
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-1 pItt ^max 

QTix) = -j= / d^Qsinm^ = ^ (-F,'"(a;)e|„ - ^T^^Wim) , (44) 
UTix) = ^ / d<pU cos = ^ (*r(a;)e|„ + Fr{x)bU , (45) 

V^TT Jo 



6.1.4. Spin-weight functions 

Functions Fp and do not form an orthogonal basis and to perform the further 
decomposition of polarization we need to use the normalized spin - weight spherical 
functions which can be defined as follows: 



^tmi^) = (-l)™+'A,-„(-x), Xj^Jx) = (-l)-+^A+„(-x) . 
These functions satisfy the self-adjoint equation: 



(47) 
(48) 



d_ 

dx 



1-^2 



A^,m + ^(^ + l)A^,^ = 0, 



and, so, they are orthogonal for a given m (see, e.g., 30, 31 and references there). 
Some of this functions can be written directly: 

A+,^ = ^^,n^(-i)^(i - xf{i - x^)^ = f:^^^z. 

Kn+l.m = A^+l,m(-l)^(l - x)\l - X^)^ [2 + (m + l)x] = F^^^ + ^^^^ , 



m(m- l)(2m+ 1)!! 



2m + 3 



(m + 2)(2m+2)!! ' ^+^'^ (m + 3)(m-l) ^'^ ' 
and using these relations functions Xi^rn with £ > m-\-l can be evaluated recursively: 



A. 



£,m 



£{£-!) 



Ci^m — 



-A^-2,n 



(53) 



' P{4£^ - 1) 
(£2 -rn2)(£2 _4) • 



dx[Fr(:r)Fr(:r) + $r(^)*r(^)] = : 



dx[Ft{x)<^f{x) + i.r(x)Fr(x)] =0. 
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This means that coefficients for the decomposition of Q and U can be found from 
(43-46) as follows: 

elm = dx[^T{x)UT{x) + Fr{x)Q^{x)\ , (57) 

e|„ = J'^ dxi^Ti^WTi^) - Fr{x)QT{x)] , (58) 

blm = dx[Fr{x)U^{x) - ^T{x)Q^{x)] , (59) 
btm = dx[-Fp{x)Ur{x) - <^Ti^)QTi^)] , 



(60) 

For numerical analysis of the polarization maps together with the temperature 
fluctuations it is convenient to use relations (27 - 38) and (57 - 60) instead of the 
recursive relations (53). Inverse problem that is the construction of the polarization 
maps from coefficients eim^^im can also be solved using relations (43-46). 

6.1.5. Spectra of the anisotropy and polarization 

Four power spectra can be introduced for the temperature, electrical and magnetic 
Stokes parameters, they are: 



m I 5 



= ^ E + i^'^mf] , (62) 

m=0 

m=0 



m=0 

When averaged over the sky, the mean square temperature anisotropy is 



iAT') = Tsf:^-^Cj, (65) 



=2 



where To is the temperature of the CMB. The mean square of polarization is 
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where 



oo 



2£+l 
An 



{El) 



'i 5 



(67) 



OO 



2£ + l 
47r 





(68) 



More details can be found in, for example, Refs. 32, 33, 34, 27. 

6.2. Incorporating of the polarization to the GLESP code 

To realize algorithms described above, special procedures had been created. We 
had developed this code in parallel in the two algorithmic languages, GNU C and 
FORTAN-77. The procedures have been designed both like subroutines and package 
utilities implementing in the GLESP pixelization scheme. 

Two commands of GLESP-pol should be used for calculation of polarization, 
namely, ^polmap' and 'polalm'. ^polmap' calculates Q and U-polarization maps by 
coefficients of E- and B- polarization modes. ''Polalm^ calculates coefficients of E- 
and B-polarization modes by Q- and U-polarization maps, respectively. E- and 
B-coefficients recorded as aim- format files, and Q- and U-maps recorded like 
GLESP maps. Both procedures can be used for temperature-spherical harmonics 
transformations too. These procedures are included to the GLESP package, version 
2.0. The programs use the fast Fourier transfom FFTW-3.2.1 Ref. 35. 

The program ' difmap' is developed for standard polarization transforms to cal- 
culate a polarization angle PA = 1/2 arctan([//(5) and an intensity / = \/Q2 + U2. 
Both values can be plotted with the GLESP drawing procedure ^f2fig\ 

6.3. Accuracy restrictions 

Both temperature and polarization maps and spectra processed by the GLESP 
package are determined in main by four parameters. First of all, it is the maximal 
number of harmonics under consideration, imax (25). The second one is the num- 
berof rings used for map presentation, Nq. By definition, we must have Nq > 2imax- 
The third one is the number of pixels for each ring of the map, N^{Oi). If we like to 
use the same angular resolution in azimuthal and polar directions than, evidently, 
we must take A^0(7r/2) ^ 2No^ and N^{Oi) decreases progressively away from the 
equator. The fourth one is the Nyquist parameter, Ny, which regulates the precisions 
achieved by calculations. 

The precision depends upon the choice A^^ > 2irnax and the Nyquist parameter, 
Ny. Thus, the orthogonality of spherical harmonics Yf^ on the sphere is achieved 
only if number of pixels for each ring is at least two times larger than number of 
harmonics we use in our analysis, N^{Oi) > 2m, < m < £ (Nyquist restriction, 
N\{Oi) > m/Ny > 2m, Ny < 0.5). Since we try to keep the same pixel area for 
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different latitude, this condition is not completely satisfied anymore because rings 
with number of pixels N^{Oi) < 2m drop out from any calculations. 

Fortunately, the influence of these restrictions on the precision achieved is mod- 
erate because the contribution of polar areas is quite small for almost all harmonics 
by the following reason: 

• Spin- weighted harmonics satisfy ±2X^771 sin^~'^{0) for m > 2. This means 
that harmonics with high m are weighted negligible everywhere apart from 
the equator, where number of pixels is enough for orthogonality. 

• Harmonics with small m < rrimin = x ^min are orthogonal. 

• Therefore, it can be expected that the harmonics with m ^ rrimin are the 
most dangerous. For such harmonics rings nearby the poles have number 
of pixels smaller than 2 x m and these parts of the sphere spoil the pic- 
ture. However, the weight of these parts of the sphere is proportional to 
siri^~'^{0) and tends to be zero for ^ ^ 0, tt. 

To accelerate the computations we remain in the calculations only larger terms 
restricted by conditions \fp{x)\ > \Fp\, \^T\ — ^ where s depends upon £ and 
required precision and, so, is determined in practice. Our tests show that the choice 
s = Si for i < 500 and e = €2 ioi £ > 500 allows to moderately decrease the 
calculation time and weakly change the precision achieved. 

We can conclude, that for a given irnax the reasonable precision can be achieved 
for the parameters : 

N0 = {2.5-3)£max. N^ = 2xN0, Nmin = 9 -11, Ny = 0.5, (69) 

and 

Sl ^ 10"^^ 62 ^ 10"^. 

Evidently, the precision achieved increases for larger Ne and AT^^, however, is 
accompanied by corresponding fast growth of the calculation time. Perhaps, for 
calculations of polarization the better results can be achieved with the repixelisation 
of the map in accordance with the net of roots of the spin - weight functions 
^fm (^^) iiistead of the Legendre polynomials (5) and corresponding choice of the 
weighting coefficients (6). 

6.4. Data format 

Developing the GLESP package for polarization calculation, we should change a 
format of data representation. We introduce two types of the format describing 
^i^m-coefficients and maps. 

In the first case, we can use the standard airn- coefficients data, which contain 
index describing number of £ and m modes corresponding to the HEALPix, real and 
imaginary parts of airn • These three parameters are described by three-fields records 
of the FITS Binary Table Ref. 36. Map data are described by the three- fields Binary 
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Table FITS format containing a vector of x^ = cos^ positions, a vector of numbers 
of pixels per each layer Nff^. , and set of temperature values in each pixel recorded by 
layers from the North Pole. All these data description formats are used separately 
for each type and polarization mode data, i.e. maps for temperature anisotropy, Q 
and U-data are contained in singles files, and aim and E and B-modes coefficients 
are stored in singles file too. 

The second type of the data representation format is similar to HEALPix one. 
In this case, data are unified in 3 extensions containing maps with anisotropy, Q- 
and U-polarization data or a^^-coefficients of temperature expansions and E- and 
B-mode, respectively. Files with coefficients in GLESP and in HEALPix have the 
same format. Each of three FITS extensions of GLESP files with maps contains 
three fields described above. 

So, using two formats of data representation a user can easily select a path of 
one's data processing including or excluding any type of polarization or anisotropy 
data. 
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